Semiclassical Solution of the Quantum Hydrodynamic Equation for Trapped 

Bose-condensed Gas in the / = Case 
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In this paper the quantum hydrodynamic equation describing the collective, low energy excitations 
of a dilute atomic Bose gas in a given trapping potential is investigated with the JWKB semiclassical 
method. In the case of spherically symmetric harmonic confining potential a good agreement is 
shown between the semiclassical and the exact energy eigenvalues as well as wave functions. It 
is also demonstrated that for larger quantum numbers the calculation of the semiclassical wave 
function is numerically more stable than the exact polynomial with large alternating coefficients. 
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I. INTRODUCTION 

The main aim of this work is to study the applica- 
bility of the JWKB semiclassical quantization method 
for the solution of the quantum hydrodynamical equa- 
tion describing the collective excitations of a trapped 
Bose-condensed gas. Although the analogous equa- 
tion for homogeneous systems has been known for long, 
the equation for trapped gases has first been written 
down and solved for special confining potential only re- 
cently by S. Stringari M (see also the review articles 
@7 §1 Rl)) after the first experimental successes of realizing 
Bose-Einstein condensation with trapped atomic vapors 
H f| Q]- This quantum hydrodynamic equation, shortly 
called 'Stringari equation', has a quite similar structure 
to the ordinary three dimensional Schrodinger equation. 
The motive of the present investigation is provided by the 
fact that similarly to the ordinary Schrodinger equation, 
the solution of the Stringari equation is also only in very 
rare, special cases available in analytic form [0, |[ |[ [Tof . 
For this reason it is important to see how the semiclassical 
methods developed for the ordinary Schrodinger equation 
can be transported to the present problem, which has a 
turning point structure of unusual type. In this work we 
demonstrate it in the simplest, also analytically treatable 
case, where the confining potential is spherically symmet- 
ric and harmonic. A comparison to the exact eigenvalues 
shows that the accuracy of the semiclassical method is 
rather good, the deviation is only a small shift in the 
energy square. 

In the first section the Stringari equation is introduced, 
rewritten in dimensionless form, separated in spherical 
polar coordinates, and the analytic results are revisited. 

In the second section we summarize briefly all the nec- 
essary information we need to know about the JWKB 
method for the present purposes, and also derive how 
the form of a second order linear differential equation 
(containing no term of first order derivative) changes un- 
der a general transformation of the independent variable. 
This result will be used in the next section to eliminate 
an undesirable term from the radial Stringari equation. 



In the third section, after applying the above men- 
tioned transformation the semiclassical solution and the 
semiclassical quantization condition are derived rigor- 
ously. 

Finally the semiclassical results are compared to the 
exact analytic ones, and some further possible applica- 
tions of the JWKB method are proposed. 



II. THE STRINGARI EQUATION 
RESULTS 



ANALYTIC 



Let us consider a dilute Bose gas trapped in a con- 
fining potential V(r). If the interparticle interaction is 
dominated by weak s-wave scattering, i.e., the (magni- 
tude) of the s-wave scattering length a is much smaller 
than the average distance between the particles, and the 
temperature T <C T c is much less than the condensation 
temperature T c , i.e., the effect of the thermal cloud on the 
condensate is negligible then the 'condensate wave func- 
tion' ^(r, t) — (iff(r,t)), which is the expectation value 
of the field operator ^(r,*), satisfies a nonlinear partial 
differential equation, the so called Gross- Pitaevskii equa- 
tion|ll],[l§|l|: 



ih 



9*(r, t) 
dt 



3 |*(r,<)| 2 )vl/(r,i), 

(la) 



where M is the mass of the particles, and g 
a constant depending on the scattering length a. 
validity conditions of the equation are 



M 



LN- 



> a 



and 



The 



(lb) 



where L is the linear size of the condensate, and N is the 
number of condensed atoms. In this paper we restrict 
our attention to the a > (i.e., g > 0) case, which means 
that the interaction is repulsive. 

In the hydrodynamic formalism two real quantities 
are introduced instead of the complex wave function of 
the condensate #(r,t) = ^/ g(r,t)e tv( - r > t \ which are the 



2 



condensate density g(r, t) = tyty* describing the magni- 
tude of the complex wave function and the velocity field 
v ( r '*) = 2§^(* V ** ~ **V*) = ^V^(r,t) related to 
the phase of the condensate wave function. (For the sake 
of simplicity the space and time arguments r, t of the 
functions are omitted, if there is no danger of confusion.) 
Taking the time derivative of the density g and using the 
time dependent Gross-Pitaevskii equation (|la|) for sub- 
stituting the time derivatives of W and \&* the following 
continuity equation is obtained: 



§ + V(,v) = 0. 



(2) 



Inserting the formula W — ^/ge' tv into the Gross- 
Pitaevskii equation (|Ta|), and separating the real part of 
it, we get 

H^ = 7 JL—&J-Q--^-(V v ) 2 -V-g Q . (3) 
dt 2My/g v 2M K ^' y ^ w 

Taking the gradient of this equation, and using the ex- 
pression v — -pVy for the velocity field, we obtain that 
the equation describing the time evolution of the velocity 
field v(r, t) is 



71/ 



<9v 
~dt 



V[V + gg- 



Mv 2 



A^) = 0. (4) 



It is worth stressing that the quantum hydrodynamic 
equations (|J) and (Q) are equivalent with the time depen- 
dent Gross-Pitaevskii equation (fla|), they do not involve 
any further approximation, and they are valid in the lin- 
ear as well as in the nonlinear regimes. 

The ground state solution of the hydrodynamic equa- 
tions is characterized by the overall zero value of the ve- 
locity field v(r, i), in which case the equation (^) reduces 
to the time independent Gross-Pitaevskii equation 

^^gJF)^[-^A + V(r)+gg(r))^), (5) 

where \x is the chemical potential of the system, and 
¥(r,t) = y/g(T)exp(-{fit). 

If the number N of the atoms in the condensate is suffi- 
ciently large, i.e., Na ^S> L (where a > is the scattering 
length of the repulsive interaction between the atoms, 
and L is the linear size of the condensate) than the ki- 
netic energy term cx -^-=Ay/g is almost everywhere neg- 
ligible with respect to the other terms in equation (Q). 
(This approximation fails only in a narrow region at the 
boundary of the condensate where g — > 0.) The hydrody- 
namic equations have the following simpler form in this 
limit: 



g + V(,v) = 0, 



„ r &v Mv 2 



0. 



(6a) 
(6b) 



The validity conditions of the above applied Thomas- 
Fermi approximation, together with the initially imposed 
conditions ( [Lb| ) are 



Na > L > N^a > 0, 



and 



T<T C . 



(6c) 



It is worth noting that this approximation has a great 
relevance, since the validity conditions (^) can be well 
satisfied in experiments Is] . 

Stringari jl| investigated the time-periodic collective 
excitational solutions of the Thomas-Fermi hydrody- 
namic equations (|^) in the linear (low energy) regime. 
The Thomas-Fermi ground state go( r ) corresponds to the 
zero velocity field vo(r) = 0, in which case the solution 



00 0) 



l (n-V(r)), ff>-V(r) >0 



otherwise 



(7) 



is obtained from the equation (^) . Linearizing the equa- 
tions @ around the ground state go and Vo , (substituting 
g{r,t) = go(r) + 5g(r,t) and v(r, t) = Vo(r) + 5v(r,t) and 
keeping the first order terms) the following equations are 
obtained for the density and velocity fluctuations 5g(r, t) 
and <Sv(r, t) 



d5g _. . , 
-^+V(g 5v) 



0. 



dS v 

M— — + gV5g = 0. 
at 



(8b) 

Taking the time derivative of ( |8a] ) and using (|8fj) to elim- 
inate (the time derivative of) the velocity fluctuation <5v 
we obtain the time dependent Stringari equation 

r d 2 5g 



M ■ 



dt 2 



V((n-V)VSg) 



(9) 



and using the Ansatz Sg(r, t) 



*£(r) for the time 



dependence of the density fluctuation Sg we arrive to the 
time independent Stringari equation 

-V((/x - F(r))V|(r)) = Mw 2 £(r), (10) 

where u> is the angular frequency of the excitation, and 
£ is defined only in the region {r|/z — V(r) > 0}. (The 
' ~ ' (tilde) sign above £ just designates that the variable 
is not dimcnsionlcss.) This equation is the starting point 
of our semiclassical investigations. 

Let us suppose that the external trap potential func- 
tion is isotropic. In this case it is natural to rewrite 
equation jic| ) in spherical polar coordinates {r, ip} 
and separate the equation using the Ansatz £(r, (p) — 
%p-Yi, m (d,(p). (As usual, r = |r| and F;, m (tf,y>) is the 
spherical harmonic function, I, m E N. Because of the 
denominator r in the Ansatz, we have L 3 |£(r)| 2 <i 3 r = 
f^l Q \x( r )\ 2 dr ■) After separation, one obtains the follow- 
ing second order ordinary differential equation for the 
radial component x{r): 



U(r)x"(r) + U'(r)x'(r) + 



-1(1 + 1) - 



X(r) = 0, (11) 
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where U(r) stands for p — V(r). 

Let us rewrite the radial equation into dimensionlcss 
form! As length unit it is reasonable to choose the 
Thomas-Fermi length ro which is defined by the equa- 
tion V(ro) = p i.e. U(ro) — 0, thus the new independent 
radial variable is p = — £ [0, 1]. For the dimensionlcss 
energy (frequency) and potential we introduce the vari- 
able '~ == ^ ' ' ' ' ' 1 " I ^Vin 4"titi n{ i r\in n{ n\ — ^i+) 



respec- 



2^tqui and the function u(p) 



tively. Using these new notations, the equation ([11]) for 
the dimensionless radial function x(p) — x( r ) — x( r op) 
has the form 



u(p)x"(p) 



2s 



u'{ P )x'{p) 

2 U (P) 



1(1 + 1) 



u'(p) 



X(j>)=0. (12) 



Further we are interested in the case when the exter- 
nal trapping potential is harmonic oscillator potential. 
This assumption makes the problem not only analyti- 
cally solvable fl], |8j [| [l(J but also has experimental rel- 
evance jli], [l5|] . For the sake of simplicity, in this paper 
we restrict our attention to isotropic harmonic oscillator 
potential V(r) = iMu^r 2 . It means that the Thomas- 



Fermi length is ro 



Mui'i 



for the dimensionless energy 



holds, and the dimensionless potential function 



becomes u(p) = 1 
takes the form 



In this case the equation ll 



x"(p) 



l — p z 
2(s 2 + l) 1(1 + 1) 



x(p) = 0, 



p e [o, l] 

(13a) 



with the normalization and continuity condition 

lim ( = 0, if I = 0, 
if I £ 0. 
(13b) 



\x(p)\ 2 d P < 



p=0 



lim 

lim xW =Q 
p\o p 



This equation is an ordinary second order linear ho- 
mogeneous differential equation with three regular sin- 
gular points at 0, 1 and oo. It is easy to see that in the 
Frobenius series expansion x(p) = P a TlmenXnP 11 of the 
solution x(p) ° m y the coefficients Xn of even indices n 
are nonzero, thus the equation ( |13| ) can be further sim- 
plified by the substitution t = fr, n(t) — x(p) °f the 
independent variable: 



«"(*) - 
( 



1 



1 - 1 
1(1 + 1 



\2t(l-t) 



4f 2 



n'(t) + 
n(t) = 0, 



t £ [0, 1] (14a) 



t=o 



dt < 



lim 2t K '(t)- K (t) 
t\0 ' 

lim^M = o, 

t\,0 v* 



0, if I = 0, 
if 1^0. 
(14b) 



This equation has again three regular singularities 
at t<°) = 0, = 1 and i (oo) = oo with indices 



„(°) - 1+1 



a 



(0) 



„(i) 



0; and 



j ± | \j2(s 2 + !) + (/ + 1/2) 2 , respect ively , what means 



that the basic solutions of the equation (14a) have power- 

(0) 

law behavior n(t) cx t" 1 - 2 at the origin t = 0, and one of 
them is locally linear (i.e., n(l) < oo, k'(1) < oo), and 
the other possesses logarithmic singularity at t = 1. Be- 
cause of the conditions (14b) the physically meaningful 



„(°) 



Hp at the origin 



solutions must have the index a\ 
and they should not have logarithmic singularity at t = 1 
Inserting the Frobenius series expansion 



K(t) =t^J2 K ° ^ 



(15a) 



n=0 



into the differential equation (14a), the recursion 

3 / 3 I s 2 \ 

(n + l)(n + I + -)nn+i = [n 2 + -7i + nl + -- — \K n 

(15b) 

is obtained for the coefficients k„, with the initial con- 
dition kq 7^ 0. The regularity condition of n( t) a t t = 1 
can be satisfied by requiring the expansion ( |l5a| ) to be 
finite, i.e., that Kk = for all k > no, what results in the 
discrete energy spectrum JlJ 



s(n, I) = \j2n 2 + 2nl + 3n + l, n,!eN. (16) 

The figures ^ and ^ show the exact radial wave function 
x(p) for a. few quantum numbers. 

It is worth noting that by changing the dependent vari- 
able according to the formula n(t) = tf* 1 w(t) the equa- 
tion can be transformed to the hypergeometric form |l6[] 

. d 2 w , , , . \ dw 
*(! - t)~nr + (c-(a + b+ I t) — - abw = 17a) 
dt dt 

with parameter values 
- 2 I 



ab 



a,b = 



e 



c = l + ~, 



I + | ± W/ 2 + / + | + 2e 2 



(17b) 



and additional physical constrains 



t l+1 /' 2 \w(t)\ 2 dt < oo, 



]hn(w(t) +2tw'(t))t— = 0, if Z = 0. 



(17c) 



limw(t)t 2 =0, 

t\,o 



if I ^ 0. 
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FIG. 1: The exact (dimensionless) radial wave function xip) f° r quantum numbers n = 0, 1, 2, 3, 4, 5 and / = 0. 




FIG. 2: The exact (dimensionless) radial wave function x(p) f° r quantum numbers n = 3 and I = 0, 1, 2, 3. 



The physically meaningful solutions [satisfying (17c)] fall 
into the so called degenerate case of the hypergeometric 
equation, in which case the solution has the form w(t) = 
t a (l — t)Pp n (t) where p n (i) is a finite polynomial of degree 
ri. The condition for this case is that either a = c — b 
or b = c — a is an integer, and it reproduces the exact 
quantization condition (filf). 



III. PRELIMINARY STEPS 



In this section first the necessary formulae of the 
JWKB method are summarized, and then, for later use, 
we briefly derive how the equation subjected to JWKB 
approximation transforms under a general change of the 
independent variable. 



A. The JWKB semiclassical method 

The JWKB method, originally developed by H. Jef- 
freys 0, G. Wentzel @, H. A. Kramers |§ and L. 
Brillouin |^0|, |TJ has still remained one of the most pow- 
erful tools for the approximative solution of wave equa- 
tions since its birth in the very early days of quantum 
mechanics [^2[ |i[ |2j|. Basically it can be used for 
obtaining a global approximation to the solution of a lin- 
ear differential equation whose highest derivative is mul- 
tiplied by a small parameter, and the method gives an 
asymptotic approximation of the real solution on a cer- 
tain interval in terms of increasing powers of this param- 
eter. 

In this paper we consider only second order homoge- 
neous linear differential equations of the following type: 



S 2 y"(x) + Q{x)y{x) = 0, 



(18) 



where S 2 is the small real parameter. In this case the 



5 



JWKB Ansatz has the form thus in the first order the semiclassical solution has the 

general form 

/ i \ 

y(x) = exp 



1 \ 

■jS s (x)\, S 5 (x) = J2 sns n{x)- (19) 

Substituting this Ansatz into the differential equa- 
tion ( |f8l) and requiring the equivalence in any orders of 
5 the following equations are obtained for the multiplier 
functions S n in the exponent 0] : 



S (x) = ±i / VQ(tjdt, (20a) 

S t {x) = ~hiQ(x), (20b) 
o / \ , • f f 5Q' 2 (t) Q"(f) \ 



(20c) 



- C x exp + Sl(x) ) + c 2 exp ( - + Sl(a; ; 



= CiQ-^exp^jT VQWd*) +C 2 g- 1 / 4 exp^-^ y/Qtfjdtj , d,C 2 e 



(21) 



r 



The JWKB solution, truncated at the iV-th order 
Vn(x) = exp (| 2n=o $ n Sn(x)) is a uniformly valid ap- 
proximation on the interval / c K in the <5 — > limit, if 
the following relations [^4| hold for the succeeding terms 
in the exponent: 

jS (x) > Si(x) > <55 2 (x) > • • • > 5 N - 1 S N {x), (22a) 



for 



0, xel. (22b) 



In concrete cases these asymptotic inequalities define 
the interval / on which the JWKB approximation can 
be applied. On the regions where the approximation 
fails (in most cases at the classical turning points, where 
Q(x) = 0), usually the multiplier function Q(x) itself 
should be approximated by a simpler function, for which 
the equation ([l8]) is exactly solvable, and then the semi- 
classical solution (of the exact equation) and the exact 
solution (of the approximate equation) should be patched 
together. This procedure will be carried out in detail for 
the radial Stringari equation, (which is an equation with 
unusual types of turning points) in the [V ;h Section. 



B. 



The effect of the change of the independent 
variable 



It is well known that the applicability of the JWKB 
method highly depends on the proper choice (transfor- 
mation) of the independent variable of the differential 



equation investigated. For this reason the formulae de- 
scribing the transformation of the differential equation 
induced by a general change of the independent variable 
are stated here. 

The starting point is the equation ([l8|). Let us trans- 
form the independent variable x in the equation to the 
new variable t by the invertible transformation 



T 



X 

x- 1 



f I 

x I 



X(t), 
T(x), 



(23a) 
(23b) 



which is defined on a suitable interval of the real line. 
Denoting the new dependent variable by y(t) = y(X(t)), 
the original equation (Rq) has the transformed form 



S 2 y"(t) + 5 2 P(t)y\t) + Q{t)y{t) 
where the multiplier functions are 
X" , x/ T" 



0. 



(24a) 



P 



X' 



■(Hx')Y 



T ,2 



oX 



(QoX)X> 2 =(^L) ,Y. 



\n{T')oX)\ 
(24b) 



The term containing the first order derivative y'(t) can 
be eliminated by applying the 



y(t) = exp ( - 

■ z ■>z=z 



P(z)dz) y(t) = -A= (25) 
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transformation of the dependent variable. After the 
transformation the equation 

5 2 y"(t) + Q(t)y(t) =0 (26a) 

is obtained, where the new coefficient function is 

Q = (QoX)X> 2 + ^ X \ xi2 iX = 

j- 2 (Q + 5i(-\u> + \tf)))oX (26b) 



with 



T" 

H= — = (lnT')', 



(26c) 



and the new unknown function y(t) is related to the orig- 
inal unknown function y(x) via the formula 



V = 



±={yoX)=(yTy)oX. 



(27) 



Thus applying two consecutive transformations to the 
differential equation (|l8|), first an arbitrary invertible 
change of the independent variable ( p3|) and then an ap- 
propriate transformation ( P5|) of the dependent variable, 
one arrives to the equation ( J26a| ) o f the same form, but 
with different coefficient function ( 26b ). This freedom 
of the choice of the independent variable enables us to 
transform the investigated equation to a more appropri- 
ate form before applying the JWKB method. 



IV. THE SEMICLASSICAL SOLUTION OF THE 
RADIAL STRINGARI EQUATION (IN CASE OF 

I = 0) 

In this section the JWKB method, summarized in the 
previous section is applied to the (dimcnsionless radial) 
Stringari equation (Jl3|) derived in Section II. 



As we did in the previous subsection [formula (25)], we 



can get rid of the term containing the first order deriva- 
tive x'(p) m equation (13a) by applying the transfor- 



r p *—% j 

mation y(p) := e ■'~=° 1 ~*' 2 z x(p) — yl P 2 x(p)- As 
result, the equation 

5 2 y"(x) + (Q(x) + 5 2 Wi(x))y{x) =0, xE [0, 1] 

(28a) 



Q(x) 



l-x 2 



Wi(x) 



1 



(l-x 2 ) 2 



(28b) 



is obtained with the normalization and continuity condi- 
tions 



1 „2 



y jx) 

l-x 2 



dx < co, 



lim ( 2MV 

x\0 \ x ) 



0, if I = 0, 



li m KM=o, 

x\0 x r ' 



(28c) 



where the abbreviations 
P 1 



2e 2 



L = 1(1 + 1) 



(28d) 



were used. (For convenience we denoted here the inde- 
pendent variable with x instead of p.) This equation is 
the starting point of the semiclassical quantization, and 
S 2 is used as small parameter. 



A comparison of th e equation ( 28a ) to the general form 
( |l8| ) shows that (28a) does not have the desired form since 
it contains a small term 8 2 Wi(x)y(x). 

First, to eliminate the Wq(x) = ( 1 _^.2)2 P ar t of this 
term, (which does not depend on the angular momentum 
quantum number I), we apply a transformation (^3| ^5|) . 
discussed in Subsection III B| , to the equation (|28[), and 
choose the transforming function T(x) to be a solution 
of the Ricatti-equation 



1 



II' (x) - -U 2 (x) = W (x) 



1 



where II = %p- 



4" (l-x 2 ) 2 ' 
= (InT')'. It is easy to see that 



(29) 



n(aO 



2.7: 



l-x 2 



T(x) = arth(x), X(t) = tanh(t) 

(30) 



is a simple solution of the nonlinear differential equation 
(^9|), and for this choice the transformed equation (|26| ) 
has the form 

S 2 y"(t) + (Q(t)+6 2 W l (t))y(t)=0, te[0,oo) (31a) 
with 

1 r ~ rr , , -L 



Q(t) 



cosh^(i) 



Wi(t) 



sinh (t) 



(31b) 



for the transformed variable y(t) — cosh(t)y( tanh(t)) . 
Furthermore, the physically relevant solutions must obey 
the supplementary conditions 



fit) 



t= o cosh (t) 



■dt < co, 



\im(m' = 

t\0 1 



We note that the correspondence 
X (tanh(<)) = y(t) 



0, if I = 0, 
if I + 0. 
(31c) 

(32) 



holds between the original variable x(p) °f the equation 
(13a) and the transformed variable y(t) of (31a). 

To get rid of the undesirable term 8 2 Wi in the follow- 
ings we restrict our attention to the I = L = case, 
where this problem does not occur, sin ce W o(t) = 0. 

In this case the differential equation (31a) has the sim- 
pler form 



S 2 y"(t) + Q(t)y(t) = 



with 



(33a) 



7 



Q{t) 



cosh 2 (i)' 



(33b) 



where the dimensionless potential function —Q(t), de- 
picted in figure ^ has a quite unusual form. There are no 
real turning points, the function — Q(t) builds a parabolic 
potential well around the origin (for < ( < 1), and it 
approaches exponentially to zero from below a s t — > oo. 
According to the formulae (pp|) in Subsection III A , the 



action variables S n (t) (and their asymptotic behavior) in 
the semiclassical solution of equation (33) are 



S (t) = ±2i (arctan(e') - - 

S 1 (t) = -ln(coBh(t))- T 
t 

~ 2 
S2(i)~±ie* 



(i»l), (34a) 

(t»l), (34b) 
(i»l), (34c) 



Applying the formula ( pl| ) and the validity criterion (22), 
(which has the form Je^k 1 in the present special case), 
the following first order semiclassical solution is obtained 
for the equation (|33|): 



(t) = A\ ^/cosh(t) sin f — (arctan(e') 

\o V 4 

if < i < -ln(tf), 



(35a) 



4 (/2 . (2 (% 
— =e ' sin 



V2 



ifl«t«-ln(i), (35b) 



where A\ S K. (In the t ^> 1 approximation the estima- 
tion f — arctanx = f°° t^-j « f 00 ^ = I was use d 
inside the argument of the sine function, which is valid 



the 



e* > 1 limit.) 



To get a good approximative solution of the equa- 
tion (|3^) for large t > 1 values it is straightforward to ap- 
proximate the coefficient function by Qo(t) = cosI ^ ^ ~ 

4e _2t if 1 <C t- With this approximation, however, the 
obtained differential equation <5 2 2/"j(t) + Ae~ 2t y\\(t) = 
is exactly solvable; the two basic solutions yn, yjj are 



y u (t) = An J I -e 
j&(t) = Ajjlo 



and 



ifl«t, 



where Jo and Yo are the zeroth order Bessel functions of 
the first and second kind, respectively [ fllij , and Au, A^ e 
R. The function yjj(t) is linearly diverging for t — > oo val- 
ues what, according to the correspondence (|32]), causes 
the divergence of x(p) at p — > 1, so the physically mean- 
ingful solution of our problem is only yu(t). 

Thus we have obtained two approximative solutions for 
the equation (B3), which are valid in two different, but 



especially for small 6 — > values exceedingly overlap- 
ping intervals, as it is shown in figure |]. To get a semi- 
classical quantization condition the functions yi(t) [equa- 
tion (HH)] and yu(t) [equation (^6|)j should be matched 
in the overlapping region 1 -C t <C — ln(<5). At first 
sight the two functions bear no readily visible resem- 
blance of form, however, using the asymptotic expansion 



cos (z — for the Bessel function, which 



is valid for \z\ 3> 1 values fllfl], the approximation 



\fl A ^ t/2cos {r 



ifl«t«- ln(8) (37) 



is obtained, which can well be compared to the for- 
mula ( |35b| ). The mat chin g condition between the two 
asymptotic formulae ( |35b|) and ([57]) is clearly the re- 
quirements An = ivM^ 1 an< ^ S ™(f (f — e l )) = 
± cos (§e~* — ?) for every t value in the interval l<f< 
— ln(<5). This means that the sum of the arguments of the 
trigonometric functions must be nix + 7r/2, wh ere n G N 
is a nonnegative integer. Using the definition (33b) of S 
the following semiclassical energy spectrum is obtained 



£ sc (n) 



2n 2 + 3n - 



1 



1 = 0, n e N, 



(38) 



which differs from the exact spectrum ( |l6| ) only by a 
constant shift of | under the square root sign. 

With the help of the formula ( [32) th e semiclassical ra- 
dial wave functions yi(t) [equation^ 35a)] andyn(i) [equa- 
tion d36|)] can be transformed back to approximate the 
original radial function x(p) [ m equation (13|)]: 



Xi(p) = w(arth(p)) =<5(f-p 2 ) *x 

'2/ [TTp n 

an i / 

V 1-p 4 

Xii(p) = Z/n(arth(p)) = 

/ 7T<5 



(-1)" 



f 



(39a) 



(39b) 



(Here the 'normalization' condition ko = y'(0) = x'(0) = 
I was used, which implies that the multiplier constants 
in the semiclassical formulae (pq) and (pfl) are A\ = 5 



(36) and 4ii = (-l)Vf .) 



To get an insight into the accuracy of the semiclas- 
sical approximation, in figures |o] and |^ we plotted the 
exact \x{p) — K (p 2 ), equation (fL5h] and the semiclassi- 
cal wave functions [xi(p)i Xu{p) equations (|39|)] as well 

as the relative errors 'Xiieh^M and xh(p)-x(p) . f or foe 

x(p) x(p) 
lowest two radial quantum numbers n = 0,1. The ex- 
act radial function was calculated with the exact energy 
eigenvalues (K&j [I = 0] and with the 'normalization' con- 
dition Kq = y\Q) = 1j while the semiclassical functions 
Xi{p) an( i Xii(p) were calculated with the semiclassical 
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FIG. 5: a) The exact and the semiclassical (dimensionless) radial wave functions [x(/°), Xi(p) an d Xn(p)]; b) The relative error 
^j,y °f the semiclassical wave functions (where Ax is Xi — X or Xn — x) f° r radial quantum number n — 0. 
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energy spectrum ( |38[) (and with the same 'normalization' 
condition). 

The figures |[ || demonstrate that even for the low- 
est energy levels the semiclassical wave function Xi(p) 
already approximates the exact function x(p) fairly well 
on almost the whole interval [0,1]. The only failure of 
Xi(p) is that it diverges at p — 1, while the exact function 
x(p) has finite value with finite slope at this point. This 
'illness' is cured by the second approximative solution 
Xn(p) which is quite different from the exact function 
on almost the whole interval [0, 1], but in the immediate 
vicinity of p = 1, where Xi{p) deviates from x(p)i gives 
a good approximation. 

We remark that the divergence of the relative error at 
the roots of the exact function x{p*) — are due to an 
unavoidable, small shift between the roots of the exact 
and semiclassical functions. 

In the end we find it worth making a final remark upon 
the numerical stability of the calculation of the exact and 
semiclassical wave function. Figure |?j shows the semiclas- 
sical and the exact wave function, xi{p) an d x{p) f° r a 
somewhat larger but still relatively small radial quantum 
number n = 13, calculated with the help of standard 
C routines with the usual floating point accuracy. It is 
clearly visible that above a certain point po ps 0.7 the 
exact polynomial totally looses its numerical stability, 
while the semiclassical function can still be easily calcu- 
lated and serves as a very good approximation almost on 
the whole interval (apart from the II nd region (0.997, 1], 
where xu(p) becomes the more appropriate approxima- 
tion). The table § shows the ( exac t values of the) co- 
efficients Kj of the polynomial (15a) obtained from the 
recursion ( |l5b| ) in the I = 0, n = 13 case. The unexpect- 
edly large values of these coefficients cause the numerical 
instability in the calculation of the exact function x(p)- 
Therefore, in certain cases it may well be reasonable to 
use semiclassical wave functions instead of the exact ones 
in numerical calculations. 



TABLE I: The coefficients k, [equation dlH)] for I 
n = 13 



and 



i 







1. 


1 


-125.6667 


2 


4674.7998 


3 


-80807.25 


4 


785626.0625 


5 


-4756608.5 


6 


19026434. 


7 


-52005588. 


8 


98657656. 


9 


-129812704. 


10 


116213280. 


11 


-67523128. 


12 


22957864. 


13 


-3466572. 



V. DISCUSSION 

In the previous sections the JWKB semiclassical ap- 
proximation has been rigorously carried out for a very 
special case of the Stringari equation ( |l0| ) (with spher- 
ically symmetric harmonic external potential V(r) and 
I = zero angular momentum quantum number). Be- 
cause of the unusual type of the turning point, the vicin- 
ity of this point needed special treatment and careful 
analysis. A comparison of the semiclassical results to the 
known exact solutions and eigenvalues shows that the ac- 
curacy of the semiclassical approximation is quite satisfy- 
ing; there is only a small (| in dimensionless units) shift 
in the energy square, and the relative error of the semi- 
classical radial wave function is already for the lowest en- 
ergy levels under a few percent. For larger quantum num- 
bers the semiclassical wave function is numerically much 
more stable than the exact polynomial, what may make 
the use of the semiclassical functions preferable than the 
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FIG. 7: The semiclassical x~i(p) an d the exact x(p) dimensionless wave functions calculated by standard C routines for 
quantum numbers I — and n = 13. 



exact ones in certain numerical calculations. 

We remark that a simpler, more practical way of 



achieving the same first order semiclassical solution (35a) 
of the radial Stringari equation (33) in the I = case 
would have been just to insert the Ansatz (|l^) right 
into the equation (28a), without bothering us about the 
5 2 Wq(x) term. This term is so small that its presence 
disturbs only the second order term S2 (x) of the JWKB 
approximation, which does not appear in the first or- 
der semiclassic al s olution (f2lj), but it does turn up in 
the condition ( |22b| ) that defines the validity region / of 
the semiclassical approximation. So it is safer anyway to 
get rid of the term S 2 Wq(x) before applying the JWKB 
method. 

Finally we note that in the works |2(| [2?], |2f| 
the same physical problem (collective excitations in 
Bose-condensed gases in spherically symmetric harmonic 
traps) has been studied with basically different semiclas- 
sical methods. The authors of these works could handle 
also the I 7^ case. They applied the Langer modifica- 
tion |2{| 1(1 + 1) — > (l + k) , and in their formula the 
shift of the energy square is 1 in dimensionless units, as 
compared to the exact result. 



As a next step, it would be useful to extend the present 
semiclassical method also for the I ^ case, and then ap- 
ply the same techniques for the solution of the Stringari 
equation in non-harmonic potentials, which are analyti- 
cally not treatable, but which have more and more ex- 
perimental relevance with the growth of the extent of the 
condensate. 
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